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Abstract 



f~i , A molecular gas system in three dimensions is numerically studied by the energy 

J^ ' conserving molecular dynamics (MD). The autocorrelation functions for the velocity 

and the force are computed and the friction coefficient is estimated. From the com- 
parison with the stochastic dynamics (SD) of a Brownian particle, it is shown that the 
force correlation function in MD is different from the delta-function force correlation 
in SD in short time scale. However, as the measurement time scale is increased further. 



"^ \ the ensemble equivalence between the microcanonical MD and the canonical SD is 



restored. We also discuss the practical implication of the result. 

Keywords: Brownian particle, molecular dynamics, stochastic dynamics, ensemble 
equivalence 



^ ' 1. Introduction 



Since Einstein published a seminal paper on the Brownian motion in 1905, it be- 
^ , came one of the most well-established subjects in statistical physics. The motion of a 

On ' Brownian particle has been studied in many works theoretically and numerically, and 

was extended later to Levy noise. Levy flights. Levy walks, continuous time random 
walks, fraction diffusion, etc. These extensions are being used to describe complex 
phenomena, e.g., anomalous diffusive behaviors llO or the diffusion limited growth 
and aggregation mechanisms [2]. For physicists, the study of Brownian motion led to a 



m 

p^ , broad class of equations of motion containing various stochastic effects. Especially, the 

Langevin equation yD is the most representative differential equation with the stochas- 
tic random forces of the white noises. Mori [4] derived the generalized Langevin equa- 
tion with non-Markovian noises, in which the memory kernel plays an important role. 
In contrast, the original Langevin equation does not have a finite memory. 

Molecular dynamics (MD) [5] is a computer simulation method to describe atoms, 
molecules, and even stars, interacting with each others in a closed system. Since the 
interaction in molecular dynamics obeys classical mechanics, the positions and mo- 
menta of particles are calculated by numerical integration of Newtonian equations of 
motion. It is important to note that the total energy of the system in MD is conserved 
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as the system moves along the trajectory in phase space. In this regard, the funda- 
mental hypothesis of equal a priori probability in statistical mechanics ensures that the 
MD simulation basically generates the microcanonical ensemble. In the computational 
point of view, the number of particles in MD is often limited and far less than any 
real system because of the limitation in computer capacity. Even with this drawback, 
MD simulation has been proven to be an excellent approximation for investigation of a 
variety of classical quantities of real materials. 

The computational limitation by the all-particle approach in the microcanonical 
MD can be overcome if one adopts a different approach. Imagine that we can concep- 
tually divide the whole system into two subsystems: One is composed of the degrees 
of freedom we like to trace, and the other is the environment system that plays the role 
of heat reservoir. In the MD approach, we need to integrate all particles' equations 
of motion. However, if it is possible to describe interaction by environmental degrees 
of freedom as stochastic random forces to the system variables, numerical integrations 
become much lighter simply due to the reduction of the number of degrees of freedom 
we need to trace. In this stochastic dynamics (SD) approach, we only need to integrate 
equations of motion for system particles, and effects from other environmental particles 
are handled as stochastic random forces that satisfy some given statistical properties. 

In the present work, we use the MD approach with total number A^ of particles, the 
volume V of the system, and the total energy E fixed (NVE ensemble in MD) |5J. In 
this MD method, it is to be noted that the time evolution of the system is deterministic 
and only depends on initial conditions within numerical accuracy. During the MD 
simulation, one is allowed to look at a single particle (call it a Brownian particle) and 
consider every other particles as composing an environment system, applying forces to 
the Brownian particle from time to time. This simple change of view allows us to make 
connection between the MD and the SD approaches, which composes the main theme 
of the present paper. In equilibrium, the energy conserving microcanonical ensemble 
is equivalent to the energy fluctuating canonical ensemble in thermodynamic limit ^ 
in most situations, with some interesting exceptions [TJ. In the present context, the 
key question to pursue in our work is in what condition the ensemble equivalence 
between the MD and the SD approaches becomes valid, in parallel to the ensemble 
equivalence between the microcanonical and the canonical ensemble in equilibrium 
statistical mechanics. In more detail, we simulate the MD with the NVE ensemble and 
observe the motion of a Brownian particle in viewpoint of the SD. All the combined 
applied forces by other particles on the Brownian particle are interpreted as the effective 
stochastic forces, and the motion of the Brownian particle is compared with that from 
the simple Langevin equation with stochastic random force. 

The present paper is organized as follows; In Section |2] we describe the details of 
our MD simulation method and quantities to be measured. The obtained results are 
shown in Section[3]in comparison with the SD approach of a Brownian particle, which 
is followed by the summary in Section |4] 

2. Simulation Methods 

In our simulations, we use A^ molecules of the identical mass m in the presence 
of the truncated Lennard-Jones interaction called the WCA (Weeks, Chandler, and 



Andersen) potential DSD, which contains only the repulsive part of the Lennard- Jones 
interaction: 

+ e, forr<rc, 

otherwise, 




12 / ^6 
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where re = 2 ' ^^cr with the interaction length scale cr and the interaction strength e. This 
representation of potential guarantees continuity of the potential and the force at r = r^, 
i.e., VwcAl=r,- = VwcAl,=r,+ = and (dVwcA/dr)\,.=r,- = (c/ywcA/£Olr=^+ = 0. We 
make equations of motion dimensionless by choosing cr, m, e, and -\Jmcr^/€ as units of 
the length, the mass, the energy, and the time, respectively, and obtain 



^-48Zf."-k)^^- 



(2) 

where 2' denotes that the sum is over only molecules (/s) satisfying r,j = |r,j(= 
rj - r,)| < Tc. Because the force term on the right-hand-side of Eq. (I2) does not depend 
on the velocity, we can use the following numerical integration scheme: v{t + At/2) - 
v(t - At 12) + a(t)At and x(t + At) - x(t) + v(t + At 12) At with the position x, the velocity 
V, the acceleration a, and the discretized time step size Af. This method is called the 
leap-frog algorithm and it is easy to see that this is the second-order algorithm although 
it runs at the same speed as in the simple first-order Euler method [[51]. 

In MD simulations, a three-dimensional cubic box of the linear size L (the vol- 
ume V - L^) with L = 24 (in unit of cr) under periodic boundary condition is used. 
The total number of particles is set lo N - 864 (the number density is thus fixed to 
p = N IV - 1/16) to make initial positions of particles fit to face-centered cubic (fee) 
structure, in order to avoid the huge amplitude of the force which might happen if the 
particles are scattered initially at random positions. The discretized time step in nu- 
merical integration is Af = 2 x 10""*, which is small enough so that further decrease 
of Af does not change results reported in the present paper The total simulation time 
is 2 X 10* Af = 400. We also verify that the total energy E of the system is conserved 
within numerical accuracy. The total momentum should also be conserved and can be 
set to zero for convenience. 

The equipartition theorem (1/2) 2j^j m,v^ = (3/2)NkBT with the Boltzmann con- 
stant kg is used to calculate the temperature T. We first generate uniform random 
velocities in [-1,1] and shift them to make the total momentum (or the velocity of the 
center of mass) zero. The initial temperature T' is then computed from the equipar- 
tition theorem, and we scale the velocities according to v, — > v,- yjT jT' to tune the 
system at the given temperature T. As time proceeds, the system approaches equilib- 
rium and the temperature is computed (call it T") again from the equipartition theorem. 
The deviation from the input temperature T is removed by the second velocity scaling 
V; -^ V; y]T IT" . After this procedure, we confirm that the temperature does not deviate 
much from the input temperature T . In our NVE ensemble MD simulations, the tem- 
perature is the control parameter and we use T - 1, 2, 3, 4, 5 and 6 in units of e/ks. It 
is then straightforward to numerically integrate the equations of motion (|2]i and all the 
presented results are obtained from the average over 100 independent runs. The key 
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Figure 1: (a) The velocity C,,(?) and (b) the force Cf(t) autocorrelation functions versus time t. The velocity 
autocorrelation Cy(t) decays exponentially in time in accord with the result from the Langevin equation for 
a Brownian particle. [Note that C,,(/) in (a) is in the log scale]. In (b), C/(f) crosses the horizontal axis at the 
cutoff time tc. The friction coefficient y,, is computed from Cy(t) via the Einstein approach, which is then 
compared with the fiiction coefficient yf from Cf{t) (see text and Table[T). 



quantities we measure during simulation is the force f, [the right-hand-side of Eq. (|2]i] 
and the velocity v,, which are then used to compute the velocity and the force auto- 
correlations defined by C,(t) = {vfHOvfHf + r)> and Cf(t) = {f^"\t')f^''\t' + t)), 
respectively, with / = \,2,- ■ -N, a - x,y,z, and (■ ■ ■ ) being the average over particles 
(0, directions {a), and time {t'), after a sufficiently long equilibration time (see Fig.lTJ. 
Before we delve into the interpretation of our MD results, we briefly review the 
Einstein approach for the stochastic Langevin equation of a Brownian particle ||3l|9|]: 
mr - -yr + i/, where m is the mass of the particle, y is the coefficient of the viscous 
friction, and // is the stochastic random force assumed to be Gaussian white noise. In 
the same dimensionless units as in our MD, the Langevin equation is written as 



r - -yr + i]. 



(3) 



with y and rj are in units of y/me/cr^ and e/o", respectively, and thus the noise correla- 
tion takes the form of 

{il^"\t')if"\t' + t)) = 2yT6{t) (4) 



for each component a - x,y, z, with T and t in units of e/kB and yjmcr- le. Following 
the Einstein approach 121 19|], it is straightforward to calculate the velocity autocorrela- 
tion function Cv{t) - (ksT/nije^'*''^'", which is written as 



C„(0 = Te-^', 



(5) 



in our dimensionless units. 



3. Results 

Fig. [Tfa) obtained from our MD simulation of Eq. (|2]l clearly shows that C,,(f = 
0) = r and that Cv{t) decays exponentially in time, which are in perfect agreement 
with the result from the Einstein approach in Eq. (|5). A simple curve fitting of Cv(t) 
to the exponential function gives us the friction coefficient y,. (we use the subscript v 



Table 1 : Friction coefficients computed from tlie velocity autocorrelation and the force autocoiTelation func- 
tions [y,. and yy in Eqs. (5) and (6), respectively] at various temperatures T. The time integration of the 
correlation function GmdU) for the MD in Eqs. (8) and ^5) is also shown together with the cutoif time t^ 
defined from zero crossing of Cf(t) in Fig. [TJb). The relaxation time scale t of GMoit) is also listed (see 
text). 
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r/ 


tc 


£ GMD(t)dt 
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1.0 


0.30(1) 


0.30(1) 


0.118 


0.97(1) 


0.048 


2.0 


0.38(1) 


0.39(1) 


0.096 


1.00(1) 


0.038 


3.0 


0.46(1) 


0.46(1) 


0.086 


0.97(1) 


0.033 


4.0 


0.51(1) 


0.52(1) 


0.078 


0.98(1) 


0.029 


5.0 


0.55(1) 


0.57(1) 


0.072 


1.00(1) 


0.026 


6.0 


0.58(1) 


0.60(1) 


0.068 


1.00(1) 


0.025 



to indicate that it is computed from the velocity correlation), which are tabulated in 
Tabled 

Alternatively, the friction coefficient can also be found from the Green-Kubo for- 
mula for the force autocorrelation function lliol lllll : 



r/ = ^ f\t(0)-m}dt (6) 

J J Jo 

in dimensionless form, where 3 in the denominator comes from the dimensionality. 
The cutoff time tc can be takenjarge enough so that the force autocorrelation function 
has arrived at plateau region lllOll . Lagar'kov and Sergeev [112ll have proposed another 
practical solution in which t^ is taken as the first zero of the force autocorrelation func- 
tion. We use the latter approach and compute the friction coefficient jf based on Eq. (|6|l 
and present results in Table [1] It is to be noted that the two different ways of comput- 
ing the friction coefficient, one from velocity autocorrelation and the other from force 
autocorrelation, give us the identical results within numerical accuracy. The agree- 
ment between y,, and yf can also be interpreted as implying the ensemble equivalence 
between the MD and the SD, since y,, is based on the expression from the Einstein 
approach for the SD, while the expression for jf is based on the MD. 

In order to check the equivalence between the MD and the SD in more detail, we 
next study the force autocorrelation function. For this, we first start from Eq. (13) and 
write 

</(0)/(0) = <[-rv'(0) + 7,(0)][-rv(f) + 77(0]> 

= r'<v(0)v(0) + </7(0)7,(0) (7) 

for one direction (we have skipped the index a for direction for brevity), and the time t 
is measured after equilibration so that the correlation is invariant under time translation. 
We define another correlation function G{t) as 

^^^^ _ <'7(0)'7(0) ^ </(0)./(0)-r'<v(0)v'(0) (g^ 

■yT jT 

In MD, we know all velocities and forces at each time, which are used to compute G{t) 
in Eq. ([8j [we call it GMoit)], combined with the friction coefficient computed above 
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Figure 2: The correlation function Gmoit) versus time t in Eq. (8) calculated for the MD simulation. Guoit) 
decays in time with finite relaxation time (see text and Table[T), which dift'ers from Gsoit) = 2(5(0- 



(Table [TJ. On the other hand, the corresponding correlation function Gsoit) for SD is 
written as Gsoit) - 25(0 from Eqs. dUi and dHJ. 

Fig. |2] shows Gufoit) at different temperatures. The relaxation time scale t for 
GMoit) can be estimated by using the method in Ref. Ii 1 31 with the time integration up 
to tc (see Eqs. (8)-(10) in Ref. llisll '). The resulting values of r at various temperatures 
are tabulated in Table[T] If the observation time scale is much larger than the relaxation 
time scale t of MD, we expect one can approximate Gmd(0 ~ Gsoit) - 26it). To 
confirm the equivalence between MD and SD in such a long-time scale, Gmoit) needs 
to satisfy 



f'GMoit)dt^ f'Gsoit)dt^2 f^6it)dt. 
Jo Jo Jo 



1, 



(9) 



where the last equality comes from the evenness of the delta function in time, and in 
the same spirit as in Eq. ^ we have used tc (see Table[T]) as the cutoff of the integration 
in Eq. (|9]l- We find that the equivalence between the MD and the SD is convincingly 
borne out as listed in Table [T] where C' GMDit)dt ^ 1 at all temperatures, as expected. 
In comparison to the MD, the use of SD has a significant benefit in practical point 
of view due to the small number of degrees of freedom to integrate. The key question 
to answer is then what is the condition for the equivalence between the MD and the 
SD. We find above that the observation time scale must be sufficiently large so that the 
correlation function Gmoit) can be approximated as a delta function as in Gsoit) to 
see the consistency between the SD and the MD: if we are only interested in long-time 
behavior, it is reasonable to use the SD instead of the MD. In other words, although the 
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Figure 3: The position of a particle in x direction versus time ; is sliown for (a) the MD and (b) the SD in 
a long-time scale of O(IO^). (c) and (d) display the velocity v(;) for the MD and the SD, respectively, in a 
shorter-time scale of 0(10). When the observation time scale is long enough, trajectories from the MD and 
the SD are qualitatively the same [see (a) and (b)] . On the other hand, in a shorter-time scale, they behave 
very differently [see (c) and (d)]. 
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Figure 4: The probability density function P(F) for the time averaged force F for the time interval fo- As ?o 
is increased [from (a) to (b)], P(F) approaches the Gaussian form, again implying the ensemble equivalence 
between the MD and the SD. 



motion of the particle in the MD in short-time scale is very different from that of the 
simple Langevin dynamics, the two cannot be distinguished in long-time scale. 

Fig.[3]displays trajectories x{t) for (a) the MD and (b) the SD at the same tempera- 
ture r = 1 in a long-time scale. In our SD simulation, we integrate the equation of mo- 
tion for a Brownian particle by using the Runge-Kutta second-order algorithm with the 
same discrete time step A = 2 x 10"'* in dimensionless time unit. As expected, the two 
trajectories look qualitatively the same. In contrast, v(f) in a shorter-time scale for (c) 
the MD and (d) the SD look quite different: for MD with the WCA potential, particles 
interact only when the distance between them is smaller than r^, and thus v(f) changes 
in time in a step-like fashion. The ensemble equivalence between the MD and the SD 
in a long-time scale displayed in Fig.[3]can also be seen in the probability density func- 
tion P{F~) for the time-averaged force F in Fig. HI where F{i) = (1/fo) f ° f{t')dt'. 
As the measurement time scale fo becomes larger, P{F) approaches the Gaussian dis- 
tribution, again revealing the ensemble equivalence between the MD and the SD in a 
long-time scale. The similar ensemble equivalence has also been discussed as the mass 
ratio between the Brownian particle and the environment particle is varied il4ll . 



4. Summary 

In summary, we have studied the gas system with the WCA potential within the 
MD approach and compared the results with the SD based on a simple Langevin equa- 
tion. The velocity and the force autocorrelation functions have been computed and a 
good agreement has been observed in the friction coefficient calculated independently 
from each correlation function. It has been revealed that as the observation time scale 
becomes much larger than the relevant relaxation time scale of correlation function, 
the ensemble equivalence between the microcanonical MD and the canonical SD ap- 
proaches is established. It is to be noted that the present study and the results are limited 
by the relatively low particle density. As the particle density is increased, the system 
will become liquid-like. In this liquid regime, the autocorrelations become more ex- 
tended in time, and the generalized Langevin formulation with memory needs to be 
used. 
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